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Abstract 

Many recent Markov chain Monte Carlo (MCMC) samplers leverage continuous 
dynamics to define a transition kernel that efficiently explores a target distribution. 
In tandem, a focus has been on devising scalable variants that subsample the data 
and use stochastic gradients in place of full-data gradients in the dynamic simu¬ 
lations. However, such stochastic gradient MCMC samplers have lagged behind 
their full-data counterparts in terms of the complexity of dynamics considered 
since proving convergence in the presence of the stochastic gradient noise is non¬ 
trivial. Even with simple dynamics, significant physical intuition is often required 
to modify the dynamical system to account for the stochastic gradient noise. In this 
paper, we provide a general recipe for constructing MCMC samplers—including 
stochastic gradient versions—based on continuous Markov processes specified via 
two matrices. We constructively prove that the framework is complete. That is, 
any continuous Markov process that provides samples from the target distribution 
can be written in our framework. We show how previous continuous-dynamic 
samplers can be trivially “reinvented” in our framework, avoiding the compli¬ 
cated sampler-specific proofs. We likewise use our recipe to straightforwardly 
propose a new state-adaptive sampler: stochastic gradient Riemann Hamiltonian 
Monte Carlo (SGRHMC). Our experiments on simulated data and a streaming 
Wikipedia analysis demonstrate that the proposed SGRHMC sampler inherits the 
benefits of Riemann HMC, with the scalability of stochastic gradient methods. 


1 Introduction 

Markov chain Monte Carlo (MCMC) has become a defacto tool for Bayesian posterior inference. 
However, these methods notoriously mix slowly in complex, high-dimensional models and scale 
poorly to large datasets. The past decades have seen a rise in MCMC methods that provide more effi¬ 
cient exploration of the posterior, such as Hamiltonian Monte Carlo (HMC) [8, 12] and its Reimann 
manifold variant [10]. This class of samplers is based on defining a potential energy function in 
terms of the target posterior distribution and then devising various continuous dynamics to explore 
the energy landscape, enabling proposals of distant states. The gain in efficiency of exploration often 
comes at the cost of a significant computational burden in large datasets. 

Recently, stochastic gradient variants of such continuous-dynamic samplers have proven quite useful 
in scaling the methods to large datasets [17, 1, 6, 2, 7]. At each iteration, these samplers use data 
subsamples—or minibatches —rather than the full dataset. Stochastic gradient Langevin dynamics 
(SGLD) [17] innovated in this area by connecting stochastic optimization with a first-order Langevin 
dynamic MCMC technique, showing that adding the “right amount” of noise to stochastic gradient 
ascent iterates leads to samples from the target posterior as the step size is annealed. Stochastic 
gradient Hamiltonian Monte Carlo (SGHMC) [6] builds on this idea, but importantly incorporates 
the efficient exploration provided by the HMC momentum term. A key insight in that paper was that 
the naive stochastic gradient variant of HMC actually leads to an incorrect stationary distribution 
(also see [4]); instead a modification to the dynamics underlying HMC is needed to account for 
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the stochastic gradient noise. Variants of both SGLD and SGHMC with further modifications to 
improve efficiency have also recently been proposed [1, 13, 7]. 

In the plethora of past MCMC methods that explicitly leverage continuous dynamics—including 
HMC, Riemann manifold HMC, and the stochastic gradient methods—the focus has been on show¬ 
ing that the intricate dynamics leave the target posterior distribution invariant. Innovating in this 
arena requires constructing novel dynamics and simultaneously ensuring that the target distribution 
is the stationary distribution. This can be quite challenging, and often requires significant physical 
and geometrical intuition [6, 13, 7]. A natural question, then, is whether there exists a general recipe 
for devising such continuous-dynamic MCMC methods that naturally lead to invariance of the target 
distribution. In this paper, we answer this question to the affirmative. Furthermore, and quite im¬ 
portantly, our proposed recipe is complete. That is, any continuous Markov process (with no jumps) 
with the desired invariant distribution can be cast within our framework, including HMC, Riemann 
manifold HMC, SGLD, SGHMC, their recent variants, and any future developments in this area. 
That is, our method provides a unifying framework of past algorithms, as well as a practical tool for 
devising new samplers and testing the correctness of proposed samplers. 

The recipe involves defining a (stochastic) system parameterized by two matrices: a positive 
semidefinite diffusion matrix, D(z), and a skew-symmetric curl matrix, Q(z), where z = (0,r) 
with 9 our model parameters of interest and r a set of auxiliary variables. The dynamics are then 
written explicitly in terms of the target stationary distribution and these two matrices. By varying 
the choices of D(z) and Q(z), we explore the space of MCMC methods that maintain the correct 
invariant distribution. We constructively prove the completeness of this framework by converting a 
general continuous Markov process into the proposed dynamic structure. 

For any given D(z), Q(z), and target distribution, we provide practical algorithms for implement¬ 
ing either full-data or minibatch-based variants of the sampler. In Sec. 3.1, we cast many previous 
continuous-dynamic samplers in our framework, finding their D(z) and Q(z). We then show how 
these existing D(z) and Q(z) building blocks can be used to devise new samplers; we leave the 
question of exploring the space of D(z) and Q(z) well-suited to the structure of the target distribu¬ 
tion as an interesting direction for future research. In Sec. 3.2 we demonstrate our ability to construct 
new and relevant samplers by proposing stochastic gradient Riemann Hamiltonian Monte Carlo , the 
existence of which was previously only speculated. We demonstrate the utility of this sampler on 
synthetic data and in a streaming Wikipedia analysis using latent Dirichlet allocation [5]. 

2 A Complete Stochastic Gradient MCMC Framework 

We start with the standard MCMC goal of drawing samples from a target distribution, which we take 
to be the posterior p(0\S) of model parameters 6 e R d given an observed dataset S. Throughout, 
we assume i.i.d. data x ^ p(x|0). We write p(6\S) oc exp(—!/(#)), with potential function 
U(6) = — X!xe<s l°gK x |0) — log p(6). Algorithms like HMC [12, 10] further augment the space 
of interest with auxiliary variables r and sample fromp(z|<S) oc exp(—iT(z)), with Hamiltonian 

H(z) = H(6,r) = U(6) + g(0,r), such that J exp(— g(Q, r))dr = constant. (1) 

Marginalizing the auxiliary variables gives us the desired distribution on 6. In this paper, we gener- 
ically consider z as the samples we seek to draw; z could represent 0 itself, or an augmented state 
space in which case we simply discard the auxiliary variables to perform the desired marginalization. 

As in HMC, the idea is to translate the task of sampling from the posterior distribution to simulating 
from a continuous dynamical system which is used to define a Markov transition kernel. That is, 
over any interval h, the differential equation defines a mapping from the state at time t to the state 
at time t + h. One can then discuss the evolution of the distribution p( z, t) under the dynamics, as 
characterized by the Fokker-Planck equation for stochastic dynamics [14] or the Liouville equation 
for deterministic dynamics [20] . This evolution can be used to analyze the invariant distribution of 
the dynamics, p s { z). When considering deterministic dynamics, as in HMC, a jump process must 
be added to ensure ergodicity. If the resulting stationary distribution is equal to the target posterior, 
then simulating from the process can be equated with drawing samples from the posterior. 

If the stationary distribution is not the target distribution, a Metropolis-Hastings (MH) correction 
can often be applied. Unfortunately, such correction steps require a costly computation on the entire 
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dataset. Even if one can compute the MH correction, if the dynamics do not nearly lead to the 
correct stationary distribution, then the rejection rate can be high even for short simulation periods 
h. Furthermore, for many stochastic gradient MCMC samplers, computing the probability of the 
reverse path is infeasible, obviating the use of MH. As such, a focus in the literature is on defining 
dynamics with the right target distribution, especially in large-data scenarios where MH corrections 
are computationally burdensome or infeasible. 

2.1 Devising SDEs with a Specified Target Stationary Distribution 

Generically, all continuous Markov processes that one might consider for sampling can be written 
as a stochastic differential equation (SDE) of the form: 

dz = f(z)d£ + v / 2D(z)dW(t), (2) 

where f(z) denotes the deterministic drift and often relates to the gradient of H( z), W(£) is a d- 
dimensional Wiener process, and D(z) is a positive semidefinite diffusion matrix. Clearly, however, 
not all choices of f(z) and D(z) yield the stationary distribution p s ( z) oc exp(— 

When D(z) = 0, as in HMC, the dynamics of Eq. (2) become deterministic. Our exposition focuses 
on SDEs, but our analysis applies to deterministic dynamics as well. In this case, our framework— 
using the Liouville equation in place of Fokker-Planck—ensures that the deterministic dynamics 
leave the target distribution invariant. For ergodicity, a jump process must be added, which is not 
considered in our recipe, but tends to be straightforward (e.g., momentum resampling in HMC). 

To devise a recipe for constructing SDEs with the correct stationary distribution, we propose writing 
f(z) directly in terms of the target distribution: 


CL r\ 

f(z) = — [D(z) + Q(z)] V-ff(z) + r(z), IMz)=^— (Dy(z) + Qy(z)). (3) 

3 = 1 3 

Here, Q(z) is a skew-symmetric curl matrix representing the deterministic traversing effects seen 
in HMC procedures. In contrast, the diffusion matrix D(z) determines the strength of the Wiener- 
process-driven diffusion. Matrices D(z) and Q(z) can be adjusted to attain faster convergence to 
the posterior distribution. A more detailed discussion on the interpretation of D(z) and Q(z) and 
the influence of specific choices of these matrices is provided in the Supplement. 

Importantly, as we show in Theorem 1, sampling the stochastic dynamics of Eq. (2) (according 
to ltd integral) with f(z) as in Eq. (3) leads to the desired posterior distribution as the stationary 
distribution: p s { z) oc exp(— H(z)). That is, for any choice of positive semidefinite D(z) and skew- 
symmetric Q(z) parameterizing f(z), we know that simulating from Eq. (2) will provide samples 
from p(0 | S ) (discarding any sampled auxiliary variables r) assuming the process is ergodic. 

Theorem 1. p s ( z) oc exp(— H(z)) is a stationary distribution of the dynamics ofEq. (2) if f(z) is 
restricted to the form ofEq. (3), with D(z) positive semidefinite and Q(z) skew-symmetric. If D(z) 
is positive definite, or if ergodicity can be shown, then the stationary distribution is unique. 


Proof. The equivalence of p s ( z) and the target p( z | S ) oc exp(— H(z)) can be shown using the 
Fokker-Planck description of the probability density evolution under the dynamics of Eq. (2) : 

d t p(z,t) = A(f.( z )p( Z;t )) + y^_A_( D . j .( z )p( z ^)). ( 4 ) 

i 1 i,j 1 3 

Eq. (4) can be further transformed into a more compact form [19, 16]: 

d t p{z, t) =V T • ( [D(z) + Q(z)] [p( z, t) ViJ(z) + Vp(z, t)] ). (5) 

We can verify thatp(z | S ) is invariant under Eq. (5) by calculating [e~ H ^\/H( z) + = 

0. If the process is ergodic, this invariant distribution is unique. The equivalence of the compact form 
was originally proved in [16]; we include a detailed proof in the Supplement for completeness. □ 
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Figure 1: The red space represents the set of all continuous Markov 
processes. A point in the black space represents a continuous 
Markov process defined by Eqs. (2)-(3) based on a specific choice of 
D(z), Q(z). By Theorem 1, each such point has stationary distribution 
p s ( z) sic p (z | S). The blue space represents all continuous Markov 
processes withp s (z) = p(z \ S). Theorem 2 states that these blue and 
black spaces are equivalent (there is no gap, and any point in the blue 
space has a corresponding D(z), Q(z) in our framework). 


2.2 Completeness of the Framework 


An important question is what portion of samplers defined by continuous Markov processes with 
the target invariant distribution can we define by iterating over all possible D(z) and Q(z)? In 
Theorem 2, we show that for any continuous Markov process with the desired stationary distribution, 
p s (z), there exists an SDE as in Eq. (2) with f (z) defined as in Eq. (3). We know from the Chapman- 
Kolmogorov equation [9] that any continuous Markov process with stationary distribution p s (z) can 
be written as in Eq. (2), which gives us the diffusion matrix D(z). Theorem 2 then constructively 
defines the curl matrix Q(z). This result implies that our recipe is complete. That is, we cover all 
possible continuous Markov process samplers in our framework. See Fig. 1. 


Theorem 2. For the SDE of Eq. (2), suppose its stationary probability density function p s ( z) 


uniquely exists, and that 


fi(z)p s (z) - £i=l ^( D ij(z)p*(z)) 


is integrable with respect to 


the Lebesgue measure, then there exists a skew-symmetric Q(z) such that Eq. (3) holds. 


The integrability condition is usually satisfied when the probability density function uniquely exists. 
A constructive proof for the existence of Q(z) is provided in the Supplement. 


2.3 A Practical Algorithm 

In practice, simulation relies on an e-discretization of the SDE, leading to a full-data update rule 

Zi +1 <- Z t -e t [(D(z t ) +Q(z t ))V-ff(z t ) +r(z t )] + Af(0, 2e t D(z t )). (6) 

Calculating the gradient of H (z) involves evaluating the gradient of U(0). For a stochastic gradient 
method, the assumption is that U(6) is too computationally intensive to compute as it relies on a sum 
over all data points (see Sec. 2). Instead, such stochastic gradient algorithms examine independently 
sampled data subsets S C S and the corresponding potential for these data: 

U{0) = -|£| y^logp(x|6»)-logp(6>); <Sc<S. (7) 

' ' xe5 


The specific form of Eq. (7) implies that U(0) is an unbiased estimator of U(0). As such, a gradient 
computed based on U(6 )—called a stochastic gradient [15] —is a noisy, but unbiased estimator 
of the full-data gradient. The key question in many of the existing stochastic gradient MCMC 
algorithms is whether the noise injected by the stochastic gradient adversely affects the stationary 
distribution of the modified dynamics (using Vf7 (0) in place of V£/(#)). One way to analyze the 
impact of the stochastic gradient is to make use of the central limit theorem and assume 

VU(0) = VU{0) + AT(0, V(0)), (8) 


resulting in a noisy Hamiltonian gradient ViJ(z) = ViJ(z) + [A/*(0, V(#)), 0] T . Simply plugging 
in Viif(z) in place of Vi7(z) in Eq. (6) results in dynamics with an additional noise term (D(z t ) + 
Q(z t )) [jV( 0, V(#)), 0] T . To counteract this, assume we have an estimate of the variance of this 

additional noise satisfying 2D(z t ) — >: 0 (i.e., positive semidefinite). With small e, this is 

always true since the stochastic gradient noise scales down faster than the added noise. Then, we 
can attempt to account for the stochastic gradient noise by simulating 


z t+1 z t — Ct 


(D(z t ) + Q(z t ))VH(z t )+r(z t ) 


+ N (0> e t(2D(z t ) — e t B t )). 


(9) 


This provides our stochastic gradient—or minibatch — variant of the sampler. In Eq. (9), the noise 
introduced by the stochastic gradient is multiplied by e t (and the compensation by e|), implying that 
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the discrepancy between these dynamics and those of Eq. (6) approaches zero as e t goes to zero. As 
such, in this infinitesimal step size limit, since Eq. (6) yields the correct invariant distribution, so 
does Eq. (9). This avoids the need for a costly or potentially intractable MH correction. However, 
having to decrease e t to zero comes at the cost of increasingly small updates. We can also use a finite, 
small step size in practice, resulting in a biased (but faster) sampler. A similar bias-speed tradeoff 
was used in [11, 3] to construct MH samplers, in addition to being used in SGLD and SGHMC. 


3 Applying the Theory to Construct Samplers 

3.1 Casting Previous MCMC Algorithms within the Proposed Framework 

We explicitly state how some recently developed MCMC methods fall within the proposed frame¬ 
work based on specific choices of D(z), Q(z) and H{ z) in Eq. (2) and (3). For the stochastic 
gradient methods, we show how our framework can be used to “reinvent” the samplers by guiding 
their construction and avoiding potential mistakes or inefficiencies caused by naive implementations. 

Hamiltonian Monte Carlo (HMC) The key ingredient in HMC [8, 12] is Hamiltonian dynamics, 
which simulate the physical motion of an object with position 0 , momentum r, and mass M on an 
frictionless surface as follows (typically, a leapfrog simulation is used instead): 

/ #£+i <— @t + 1 r t 

1 r t +i <-r t -e t V£/(0 t ). 1 ' 

Eq. (10) is a special case of the proposed framework with z = (0,r),H ( 6 , r) = U(0) + \r T M~ x r, 

Q(M=(J - I )andD(^r) = 0. 


Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) As discussed in [6], simply replac 
ing Vf7(0) by the stochastic gradient VC/ (6) in Eq. (10) results in the following updates: 

Ot+1 <- Ot + etMjVt 
n+ 1 <-r t - e t VU(0 t ) 


Naive : 


r t - e t VE/(0 t ) + A/(O,e t 2 V(0 t )), 


(ID 


where the « arises from the approximation of Eq. (8). Careful study shows that Eq. (11) cannot be 
rewritten into our proposed framework, which hints that such a naive stochastic gradient version of 
HMC is not correct. Interestingly, the authors of [6] proved that this naive version indeed does not 
have the correct stationary distribution. In our framework, we see that the noise term A/(0, 2e t D(z)) 
is paired with a D(z)ViT(z) term, hinting that such a term should be added to Eq. (11). Here, 

D(0, r) = ^ o eV(0) means we nee d t0 a dd E)(z)Vi7(z) = eV(0)V r i7(0, r) = 

eV(0)M _1 r. Interestingly, this is the correction strategy proposed in [6], but through a physical 
interpretation of the dynamics. In particular, the term eV(0)M _1 r (or, generically, CM -1 r where 
C y eV(0)) has an interpretation as friction and leads to second order Langevin dynamics: 

f 0 t+ i «- 6 t + 

l r m <r- r t - e t VU{O t ) - etCM-Vt +A/(0,e t (2C - e t B t )). 1 J 


Here, B t is an estimate of V(0t). This method now fits into our framework with H ( 0 , r) and Q(0, r) 
as in HMC, but with D(0, r) = ( ? r* ) • This example shows how our theory can be used to 


identify invalid samplers and provide guidance on how to effortlessly correct the mistakes; this is 
crucial when physical intuition is not available. Once the proposed sampler is cast in our framework 
with a specific D(z) and Q(z), there is no need for sampler-specific proofs, such as those of [6]. 


Stochastic Gradient Langevin Dynamics (SGLD) SGLD [17] proposes to use the following first 
order (no momentum) Langevin dynamics to generate samples 

Ot+i <^0 t - e t DV£/(0 t ) + A/(0,2e t D). (13) 

This algorithm corresponds to taking z = 6 with H(0) = U(6), D(0) = D, Q(0) = 0, and B t = 0. 
As motivated by Eq. (9) of our framework, the variance of the stochastic gradient can be subtracted 
from the sampler injected noise to make the finite stepsize simulation more accurate. This variant of 
SGLD leads to the stochastic gradient Fisher scoring algorithm [1]. 
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Stochastic Gradient Riemannian Langevin Dynamics (SGRLD) SGLD can be generalized to 
use an adaptive diffusion matrix D(0). Specifically, it is interesting to take D(0) = G _1 (6 > ), where 
G(0) is the Fisher information metric. The sampler dynamics are given by 

Ot +1 <r-9t - e t [G(0 t y l VU(0 t ) +T(0 t )\ +U(0,2e t G(e t )~ 1 ). (14) 


Taking D(0) = G(0) \ Q(0) = 0, and = 0, this SGRLD [13] method falls into our frame¬ 


work with correction term T^ (0) = Y 


9D ij(0) 

80, 


. It is interesting to note that in earlier literature [10], 


d 

Ti(6) was taken to be 2 |G(0)| -1 / 2 Y (G -” 1 (6^)| G (^)| 1//2 ). More recently, it was found that 

j vOj 

this correction term corresponds to the distribution function with respect to a non-Lebesgue mea¬ 
sure [18]; for the Lebesgue measure, the revised Ti(0) was as determined by our framework [18]. 
Again, we have an example of our theory providing guidance in devising correct samplers. 


Stochastic Gradient Nose-Hoover Thermostat (SGNHT) Finally, the SGNHT [7] method in¬ 
corporates ideas from thermodynamics to further increase adaptivity by augmenting the SGHMC 
system with an additional scalar auxiliary variable, £. The algorithm uses the following dynamics: 

@t+ i Ot + e tG ^ 

7^+i <— G — e t \/U(O t ) — c t £ t r t + A/"(0, e t (2A — e^B^)) 

f 1 \ ^ 

6+i 6 + e t r t ~ 1J • 


We can take z = (0, r, £), 77(0, r, £) = U(6) + ^r T r + — A) 2 , D($, r, = f o a i o 

A Act \ o o o 

( 0 -I 0 \ 

andQ(6>,r, 6= 1 o r/d to place these dynamics within our framework. 

y 0 -r T /d 0 J 


Summary In our framework, SGLD and SGRLD take Q(z) = 0 and instead stress the design of 
the diffusion matrix D(z), with SGLD using a constant D(z) and SGRLD an adaptive, 0-dependent 
diffusion matrix to better account for the geometry of the space being explored. On the other hand, 
HMC takes D(z) = 0 and focuses on the curl matrix Q(z). SGHMC combines SGLD with HMC 
through non-zero D(0) and Q(0) matrices. SGNHT then extends SGHMC by taking Q(z) to be 
state dependent. The relationships between these methods are depicted in the Supplement, which 
likewise contains a discussion of the tradeoffs between these two matrices. In short, D(z) can 
guide escaping from local modes while Q(z) can enable rapid traversing of low-probability regions, 
especially when state adaptation is incorporated. We readily see that most of the product space 
D(z) x Q(z), defining the space of all possible samplers, has yet to be filled. 

3.2 Stochastic Gradient Riemann Hamiltonian Monte Carlo 


In Sec. 3.1, we have shown how our framework unifies existing samplers. In this section, we now use 
our framework to guide the development of a new sampler. While SGHMC [6] inherits the momen¬ 
tum term of HMC, making it easier to traverse the space of parameters, the underlying geometry of 
the target distribution is still not utilized. Such information can usually be represented by the Fisher 
information metric [10], denoted as G(0), which can be used to precondition the dynamics. For our 
proposed system, we consider 77(0, r) = U (0) + \r T r , as in HMC/SGHMC methods, and modify 
the D(0, r) and Q(0, r) of SGHMC to account for the geometry as follows: 


D(0,r) 


0 0 \ 

0 G(0)- 1 J ; 


Q(0,r) = 


0 

G(0)- 1/2 


-G(0)- 1/2 

0 


We refer to this algorithm as stochastic gradient Riemann Hamiltonian Monte Carlo (SGRHMC). 
Our theory holds for any positive definite G(0), yielding a generalized SGRHMC (gSGRHMC) 
algorithm, which can be helpful when the Fisher information metric is hard to compute. 


A naive implementation of a state-dependent SGHMC algorithm might simply (i) precondition the 
HMC update, (ii) replace VC/(0) by VC/(0), and (iii) add a state-dependent friction term on the 
order of the diffusion matrix to counterbalance the noise as in SGHMC, resulting in: 

i\t • f 0t+1 0t + e*G(0t) 1//2 ?+ ^ tl bi 

aiVei \ r t+ i<- n-etG(0O“ 1/2 V 0 L(0t)-etG(0O“ 1 n+AT(O,et(2G(0O“ 1 -etB t )). 1 } 


6 




Algorithm 1: Generalized Stochastic Gradient Riemann Hamiltonian Monte Carlo 

initialize (0 q ? A)) 
for t = 0,1 , 2 do 

optionally, periodically resample momentum r as ~ A/"(0,1) 

0*+i F- + etG(6 > t) _1 / 2 rt, St F- et(2G(6 > t) _1 — e t B t ) 

r t+1 <- r t - e t G(0 t )- 1 /2 V() C/(0 t ) + e t V e (G(0 t )“ 1/2 ) - e t G(0 t )- 1 r t + Af{o, E t ) 




Figure 2: Left: For two simulated ID distributions defined by U{6) — 0 2 / 2 {one peak) and U {0) = 9 4 — 2 0 2 
{two peaks), we compare the KL divergence of methods: , SGHMC, the naive SGRHMC of Eq. (16), and 

the gSGRHMC of Eq. (17) relative to the true distribution in each scenario (left and right bars labeled by 1 and 
2). Right: For a correlated 2D distribution with U{Q i, # 2 ) = 0i/lO + (4 • {62 + 1.2) — 0 2 ) 2 / 2, we see that 
our gSGRHMC most rapidly explores the space relative to SGHMC and . Contour plots of the distribution 
along with paths of the first 10 sampled points are shown for each method. 

However, as we show in Sec. 4.1, samples from these dynamics do not converge to the desired 
distribution. Indeed, this system cannot be written within our framework. Instead, we can simply 
follow our framework and, as indicated by Eq. (9), consider the following update rule: 

f Ot+i 0t + etG(0t) 1/2 r t 

| rt +1 <- n - e t [G{e)-^VeUidt) + Ve (g (6t)~ 1/2 ) ~ G(0 t ) _1 n] + V(0, e t (2G(0 t ) _1 - e*B t )), 

(17) 

which includes a correction term (G(#) -1 / 2 ), with i-th component (G(6 > ) -1 / 2 )... The 

practical implementation of gSGRHMC is outlined in Algorithm 1. 

4 Experiments 

In Sec. 4.1, we show that gSGRHMC can excel at rapidly exploring distributions with complex 
landscapes. We then apply SGRHMC to sampling in a latent Dirichlet allocation (LDA) model on 
a large Wikipedia dataset in Sec. 4.2. The Supplement contains details on the specific samplers 
considered and the parameter settings used in these experiments. 

4.1 Synthetic Experiments 

In this section we aim to empirically (i) validate the correctness of our recipe and (ii) assess the 
effectiveness of gSGRHMC. In Fig. 2(left), we consider two univariate distributions (shown in the 
Supplement) and compare SGLD, SGHMC, the naive state-adaptive SGHMC of Eq. (16), and our 
proposed gSGRHMC of Eq. (17). See the Supplement for the form of G(0). As expected, the naive 
implementation does not converge to the target distribution. In contrast, the gSGRHMC algorithm 
obtained via our recipe indeed has the correct invariant distribution and efficiently explores the dis¬ 
tributions. In the second experiment, we sample a bivariate distribution with strong correlation. The 
results are shown in Fig. 2(right). The comparison between SGLD, SGHMC, and our gSGRHMC 
method shows that both a state-dependent preconditioner and Hamiltonian dynamics help to make 
the sampler more efficient than either element on its own. 


7 



















Parameter 0 
Prior p(0) 


Method 

SGHMC 

SGRLD 

SGRHMC 



Figure 3: Upper Left: Expanded mean parameterization of the LDA model. Lower Left: Average runtime per 
100 Wikipedia entries for all methods. Right: Perplexity versus number of Wikipedia entries processed. 


4.2 Online Latent Dirichlet Allocation 

We also applied SGRHMC (with G(0) = diag(#) -1 , the Fisher information metric) to an online 
latent Dirichlet allocation (LDA) [5] analysis of topics present in Wikipedia entries. In LDA, each 
topic is associated with a distribution over words, with /3k w the probability of word w under topic k. 
Each document is comprised of a mixture of topics, with tt^ the probability of topic k in document 
d. Documents are generated by first selecting a topic ~ i for the jth word and then drawing 
the specific word from the topic as ~ /3 (d) . Typically, and f3k are given Dirichlet priors. 

The goal of our analysis here is inference of the corpus-wide topic distributions /3k . Since the 
Wikipedia dataset is large and continually growing with new articles, it is not practical to carry out 
this task over the whole dataset. Instead, we scrape the corpus from Wikipedia in a streaming man¬ 
ner and sample parameters based on minibatches of data. Following the approach in [13], we first 
analytically marginalize the document distributions and, to resolve the boundary issue posed by 
the Dirichlet posterior of /?& defined on the probability simplex, use an expanded mean parameter¬ 
ization shown in Figure 3(upper left). Under this parameterization, we then compute V logp(0|x) 
and, in our implementation, use boundary reflection to ensure the positivity of parameters Okw • The 
necessary expectation over word-specific topic indicators zf^ is approximated using Gibbs sampling 
separately on each document, as in [13]. The Supplement contains further details. 

For all the methods, we report results of three random runs. When sampling distributions with 
mass concentrated over small regions, as in this application, it is important to incorporate geometric 
information via a Riemannian sampler [13]. The results in Fig. 3 (right) indeed demonstrate the im¬ 
portance of Riemannian variants of the stochastic gradient samplers. However, there also appears to 
be some benefits gained from the incorporation of the HMC term for both the Riemmannian and non- 
Reimannian samplers. The average runtime for the different methods are similar (see Fig. 3 (lower 
left)) since the main computational bottleneck is the gradient evaluation. Overall, this application 
serves as an important example of where our newly proposed sampler can have impact. 


5 Conclusion 


We presented a general recipe for devising MCMC samplers based on continuous Markov processes. 
Our framework constructs an SDE specified by two matrices, a positive semidefinite D(z) and a 
skew-symmetric Q(z). We prove that for any D(z) and Q(z), we can devise a continuous Markov 
process with a specified stationary distribution. We also prove that for any continuous Markov pro¬ 
cess with the target stationary distribution, there exists a D(z) and Q(z) that cast the process in our 
framework. Our recipe is particularly useful in the more challenging case of devising stochastic gra¬ 
dient MCMC samplers. We demonstrate the utility of our recipe in “reinventing” previous stochastic 
gradient MCMC samplers, and in proposing our SGRHMC method. The efficiency and scalability 
of the SGRHMC method was shown on simulated data and a streaming Wikipedia analysis. 
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Supplementary Material: A Complete Recipe for Stochastic Gradient MCMC 
A Proof of Stationary Distribution 

In this section, we provide a proof for Theorem 1. To prove the theorem, we first show when f(z) 
satisfies Eq. (3) in the main paper, then the following Fokker-Planck equation of the dynamics: 

d t p(z,t) = - ^ ^-(fj(z)p(z,f)) + ( D fi( Z MM)), 

i ^ i,j 1 3 

is equivalent to the following compact form [19, 16]: 

d t p{ z, t ) =V T • ( [D(z) + Q(z)] [p( z, f)ViJ(z) + Vp(z, t)] ). (S.l) 


Proof. The proof is a re-writing of Eq. (S.l): 

d t p( z, t) =V T ■ ( [D(z) + Q(z)] [p(z, t)VH( z) + Vp(z, f)] ) 


= E ^ E I D *i ( z ) + Qij( z ) 

* = 1 1 J 




5 


+ r + E^[ D f?( z ) + c M z )]E 


<9z 


<9z ? - 


<9z, 


We can further decompose the second term as follows 


E ^T D ^ z )E -E - E 7£T ] P'M 


dzj 


13 


dzi dzj 


dz i 


E 


E E ST?< z - •) = - E ST { ') 


dzj 


dz i 


E az Qii(z) 


9 9 

The second equality follows because JE ? . —— —— [Q^ (z)p(z, £)] = 0 due to anti-symmetry of Q. 


JlJ dzi dzn 


Putting these back into the formula, we get 




E [ D 'j( z > + Qij( z )l - E-^r D ' :z;i + #-%«] 


p(z,t) 




13 


= ~ E^:( f ‘( z )p( z >*)) + E ^r( D *i( z )p( z ’*))- 




□ 


We can then verify that p(z | 5) oc e ^( z ) is invariant under the compact form by calculating 


[p( z, t)VH(z) + Vp(z, £)] oc [e“^ (z) Vi^(z) + Ve“^ (z) 


= 0 . 


This completes the proof of Theorem 1 . 

The above proof follows directly from [16] and is provided here for readers’ convenience. 


-p(z,£). 
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B Proof of Completeness 

In this section, we provide a constructive proof for Theorem 2, the existence of Q(z). 

The proof is first outlined as follows: 

• We first rewrite Eq. (3) in the main paper and notice that finding matrix Q(z) is equivalent 
to finding the matrix Q(z)p s (z) such that 

53 (Qy( z )P S ( z )) = fi(z)p s (z) - J^Dy^tz)), 

3 3 3 3 

where the right hand side is a divergence-free vector. 

• We transform the above equation and its constraint into the frequency domain and obtain a 
set of linear equations. 

• Then we construct a solution to the linear equations and use inverse Fourier transform to 
obtain Q(z). 

The complete procedure is: 


Proof. Multiplying p s (z) on both sides of Eq. (3) in the main paper, and noting that: 

p s (z) = exp (-H (z)), 

we arrive at: 

d 


The equation for Q^-(z) can now be written as: 
d 


fj(z)p s (z) = £^r(( D «( z ) + Qii( z ))p s ( z )j- 


53 ^r(Q*i( z )p s ( z )) = f *( z )p s ( z ) - 53^( D «( z )p'( z ))- 


(S.2) 


(S.3) 


(S.4) 


Recall that the Fokker-Planck equation for the stochastic process, Eq. (2), is: 

dP fo - = -V T • (t(z)p(z,t)) + V 2 : (D(z)p(z,t)) 


i { 3 


(S.5) 


We can immediately observe that the right hand side of Eq. (S.4) has a divergenceless property by 
substituting the stationary probability density function p s (z) into Eq. (S.5): 




d 


= 0 . 


(S.6) 


The nice forms of Eqs. (S.4) and (S.6) imply that the questions can be transformed into a lin¬ 
ear algebra problem once we apply a Fourier transform to them. Denote the Fourier transform of 

Q(z)p s (z) as Q(k); and Fourier transform of f*(z)p s (z) — (d^- (z)p s (z)) as F^(k), where 

k = (ki, • • • , k n ) T is the set of the spectral variables. That is: 

Qij ( k ) = [ Qii( z )p s ( z )« 

JV 

d 


)e~ 2xi k Z dz; 


Fi(k) = l ( W(z) - 53 {°ij( z )P S ( z 


— 2-7ri k z 


dz. 
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d 

Then, — (Q ij (z)p s (z)j is 
alent form in Fourier space: 


transformed to 2tti Q ZJ k 7 , and Eq. (S.4) becomes the following equiv- 


2tti Qk = F 
k T F = 0. 


(S.7) 


Hence, it is clear that matrix Q must be a skew-symmetric projection matrix from the span of k 
to the span of F, where k and F are always orthogonal to each other. We thereby construct Q as 
combination of two rank 1 projection matrices: 


A , i Fk T , . wl kF T 

Q = ( 27rl ) V7Ti7 “ ( 27n ) 


(S.8) 


k T k v ' k T k 

We arrive at the final result that matrix Q(z) is equal to p s (z)~ 1 times the inverse Fourier transform 
of Q(k): 


Q, 


i (z) = p s (z) 1 [ 

Jx 


kjFj(k) - kjFj (k) 2 W iEk,x, 
v (27ri)-Ek? 6 


dk. 


(S.9) 


Thus, if |^(z)p s (z) — belongs to the space of L 1 , then any continuous 

time Markov process, Eq. (2), can be turned into this new formulation. □ 

Remark 1. Entries in the skew-symmetric projector Q^-(z) constructed here are real. 

k ? ; 


Denote a ‘f = kf, then the inverse Fourier transform of 


l^i 


M-Ek? 


along the partial variable 


k i is equal to: 


gi(z) = z<] 


where H[x\ is the Heaviside function. Because gi(z) is an even function in ki, l i, its total inverse 
Fourier transform is real. 


Therefore, the inverse Fourier transform of 


kjgj(k) 

M-Ek? 

I 


is the convolution of two real functions. 


C 2-D Case as a Simple Intuitive Example of the Construction 

For 2-dimensional systems, we have: 

kiFi(k)+k 2 F 2 (k) =0, (S.10) 

and hence Eq. (S.9) has a simple form: 

Q2l(zi,Z2) = -Ql2(zi,Z 2 ) 

= P s (zi,z 2 ) _1 [ f 1 (z 1 ,s)p 8 (z 1 ,s)ds 

OZ2 Q ^ 

-J ^(Dii(zi,s)p s (zi,s)) - — ^Di 2 (zi,s)p s (zi,s))ds 
f Zl 

= ~Ps( Zl,Z 2 ) _1 / f 2 (s,z 2 )p s (s,z 2 )ds 

j «; 

+ J o ^(d 21 (s,z 2 ) P '( S ,z 2 )) + A(D 22 ( s ,z 2 )^( S ,z 2 ))d S . CS.11) 
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D Previous MCMC Algorithms in the Form of Continuous Markov 
Processes as Elements in the Current Recipe 

This section parallels that of Sec. 3.1 of the main paper, but in terms of the continuous dynamics 
underlying the samplers. This allows us to rapidly draw connections with our SDE framework of 
Sec. 2.1. Fig. S.4 provides a cartoon visualization of the portion of the product space D(z) x Q(z) 
already covered by past methods, after casting these methods in our framework below. Our proposed 
gSGRHMC method covers a portion of this space previously not explored. 


Q(z) 



Figure S.4: Cartoon of how previous methods explore the space of possible D(z) and Q(z) matrices, 
along with our proposed gSGRHMC method of Sec. 3.2. 


Hamiltonian Monte Carlo (HMC) The continuous dynamics underlying Eq. (10) in the main 
paper are 


/ dO = M~ 1 rdt 
\ dr = -W(<9)dt. 


(S.12) 


Again, we see Eq. (S.12) is a special case of our proposed framework with z = ( 6 , r), H(6 , r) = 
U(0) + \r T M~ l r, Q(0, r) = f J J and D(0, r) = 0. 


Stochastic Gradient Hamiltonian Monte Carlo (SGHMC) As described in [6], replacing 
VU(6) by the stochastic gradient VU{6) in the e-discretized HMC system of Eq. (10) (resulting 
in Eq. (11)) has a continuous-time representation as: 


Naive : 


d 0 = M~ 1 rdt 

dr = - W(<9)d£ + v^)dW(t) « -VU(0)dt. 


(S.13) 


Analogously to Sec. 3.1, these dynamics do not fit into our framework. Instead, in our framework 
we see that the noise term ^2D(z)dW(t) is paired with a D(z)ViT(z) term, hinting that such a 

term must be added to the dynamics of Eq. (S.13). Here, D(0, r) = ^ q ^y which means 

we need to add a term of the form D(z)ViT(z) = eVV r i7($, r) = eVM _1 r. Interestingly, this 
is the correction strategy proposed in [6], but through a physical interpretation of the dynamics. In 
particular, the term eVM -1 r (or, generically, CM -1 r) has an interpretation as a friction term and 
leads to second order Langevin dynamics: 


d 0 = M~ 1 rdt 

dr = —VU(6)dt - CM~ 1 rdt + ^ 2 C - eV((9)dW(f) + ^/eV{6)dW{t). 


This method now fits into our framework with H(6,r) and Q(0, r) as in HMC, but here with 
D(0, r )=( ° ° 
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Stochastic Gradient Langevin Dynamics (SGLD) SGLD [17] proposes to use the following first 
order (no momentum) Langevin dynamics to generate samples 

d9 = -T>VU(0)dt + V2DdW(t). (S.15) 

This algorithm corresponds to taking z = 6 with H(0) = U(0 ), D(0) = D, Q (6) = 0. As in the 
case of SGHMC, the variance of the stochastic gradient can be subtracted from the sampler injected 
noise V2DW (t) to make the finite stepsize simulation more accurate. This variant of SGLD leads 
to the stochastic gradient Fisher scoring algorithm [1]. 

Stochastic Gradient Riemannian Langevin Dynamics (SGRLD) SGLD can be generalized to 
use an adaptive diffusion matrix D($). Specifically, it is interesting to take D($) = G _1 (6 > ), where 
G(0) is the Fisher information metric. The sampler dynamics is given by 

d0 = -G~ 1 (9)VU(9)dt + T(9) + v /2G irT (0)dW(i). (S.16) 

Taking D(0) = G _1 (#) and Q(0) = 0, the SGRLD method falls into our current framework with 

9D • • (0) 

the correction term Ti(0) = —ttx-• 

j Wj 


Stochastic Gradient Nose-Hoover Thermostat (SGNHT) Finally, the continuous dynamics un¬ 
derlying the SGNHT [7] algorithm in Sec. 3.1 are 

d 0 rdt 


dr = —WU(6)dt 

= (y rTr - ] 


d t + y/2A dW (t) 
rdt. 


(S.17) 


Again, we see we can take z = 
/ o o o \ 

( o a- i o J, and Q(0,r,g) 

work. 


(9, r, 0. H(6, r, £) = U{9) + l -r T r +±(£-A) 2 , D (9, r, £) = 

/ 0 -I o \ 

= | 1 o r /d J to place these dynamics within our frame- 


E Discussion of Choice of D and Q 

A lot of choices of D(z) and Q(z) could potentially result in faster convergence of the samplers than 
those previously explored. For example, D(z) determines how much noise is introduced. Hence, an 
adaptive diffusion matrix D(z) can facilitate a faster escape from a local mode if | |D(z) 11 is larger in 
regions of low probability, and can increase accuracy near the global mode if | |D(z) 11 is smaller in 
regions of high probability. Motivated by the fact that a majority of the parameter space is covered 
by low probability mass regions where less accuracy is often needed, one might want to traverse 
these regions quickly. As such, an adaptive curl matrix Q(z) with 2-norm growing with the level 
set of the distribution can facilitate a more efficient sampler. We explore an example of this in the 
gSGRHMC algorithm of the synthetic experiments (see Supp. F.l). 


F Parameter Settings in Synthetic and Online Latent Dirichlet Allocation 
Experiments 

F.l Synthetic Experiments 


In the synthetic experiment using gSGRHMC, we specifically consider G(0) 1 = Dy \ U{6) + Cj. 
The constant C ensures that U(6) + C is positive in most cases so that the fluctuation is indeed 
smaller when the probability density function is higher. Note that we define G(0) in terms of U ( 0 ) 
to avoid a costly full-data computation. We choose D — 1.5 and C = 0.5 in the experiments. The 
design of G is motivated by the discussion in Supp. E, taking Q(0) to have 2-norm growing with 
the level sets of the potential function can lead to faster exploration of the posterior. 

Comparison of SGLD, SGHMC, the naive implementation of SGRHMC (Eq. (16)), and the gS¬ 
GRHMC methods is shown in Fig. S.5, indicating the incorrectness of the naive SGRHMC. 
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Figure S.5: For two simulated ID distributions (black) defined by U (6) = 6 2 /2 (left) and U(6) = 
0 4 - 2 0 2 ( right ), comparison of , SGHMC, the naive SGRHMC of Eq. (16), and the gSGRHMC 
of Eq. (17) in the main paper. 


F.2 Online Latent Dirichlet Allocation Experiment 


In the online latent Dirichlet allocation (LDA) experiment, we used minibatches of 50 documents 
and Km 50 topics. Similar to [13], the stochastic gradient of the log posterior of the parameter 6 
on a minibatch S is calculated as 


dlogp(fl|x,a,7) 


a — 1 


7 kw 


- i 


M 'Vf 

I £Ti Z_^ ^z( rf )|x( d ),6>,7 

I I des 


dkw 

@kw 


™dk- 

Ok- 


(S.18) 


where a is the hyper-parameter for the Gamma prior of per-topic word distributions, and 7 for the 
per-document topic distributions. Here, n^kw is the count of how many times word w is assigned to 

topic k in document d (via = k for xj = w). The • notation indicates n^k- — ndkw To 
calculate the expectation of the latent topic assignment counts ridkw , Gibbs sampling is used on the 
topic assignments in each document separately, using the conditional distributions 


p( 7° = fc|x (d) ,0, 7 ) = 


( 7 +»*.) 


6 kx <<> 


E fc (■ 


'r + n dk.)o kx f 


(S.19) 


where \j represents a count excluding the topic assignment variable being updated. See [13] 
for further details. 


We follow the experimental settings in [13] for Riemmanian samplers (SGRLD and SGRHMC), 
taking the hyper-parameters of Dirichlet priors to be 7 = 0.01 and a = 0.0001. Since the non- 
Riemmanian samplers (SGLD and SGHMC) do not handle distributions with mass concentrated 
over small regions as well as the Riemmanian samplers, we found 7 = 0.1 and a = 0.01 to 
be optimal hyper-parameters for them and use these instead for SGLD and SGHMC. In doing so, 
we are modifying the posterior being sampled, but wished to provide as good of performance as 
possible for these baseline methods for a fair comparison. For the SGRLD method, we keep the 


stepsize schedule of e t 



—c 

and corresponding optimal parameters a, 6 , c used in 


the experiment of [13]. For the other methods, we use a constant stepsize because it was easier to 
tune. (A constant stepsize for SGRLD performed worse than the schedule described above, so again 
we are trying to be as fair to baseline methods as possible when using non-constant stepsize for 
SGRLD.) A grid search is performed to find e t = 0.02 for the SGRHMC method; e t = 0.01, D — I 
(corresponding to Eq. (13) in the main paper) for the SGLD method; and e t = 0.1, C = M — I 
(corresponding to Eq. (12) in the main paper) for the SGHMC method. 


For a randomly selected subset of topics, in Table S.l we show the top seven most heavily weighted 
words in the topic learned with the SGRHMC sampler. 
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“ENGINES” 

speed 

product 

introduced 

designs 

fuel 

quality 

“ROYAL” 

britain 

queen 

sir 

earl 

died 

house 

“ARMY” 

commander 

forces 

war 

general 

military 

colonel 

“STUDY” 

analysis 

space 

program 

user 

research 

developed 

“PARTY” 

act 

office 

judge 

justice 

legal 

vote 

“DESIGN” 

size 

glass 

device 

memory 

engine 

cost 

“PUBLIC” 

report 

health 

community 

industry 

conference 

congress 

“CHURCH” 

prayers 

communion 

religious 

faith 

historical 

doctrine 

“COMPANY” 

design 

production 

produced 

management 

market 

primary 

‘PRESIDENT” 

national 

minister 

trial 

states 

policy 

council 

“SCORE” 

goals 

team 

club 

league 

clubs 

years 


Table S. 1: The top seven most heavily weighted words (columns) associated with each of a randomly 
selected set of 11 topics (rows) learned with the SGRHMC sampler from 10,000 documents (about 
0.3% of the articles in Wikipedia). The capitalized words in the first column represent the most 
heavily weighted word in each topic, and are used as the topic labels. 
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